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Abstract : We consider a distributed system with persistent memory of a 
type which is often encountered in viscoelasticity or in the study of diffusion 
processes with memory. The relaxation kernel, i.e. the kernel of the memory 
term, is scarcely known from hrst principles, and it has to be inferred from 
experiments taken on samples of the material. We prove that two boundary 
measures give a linear Volterra integral equation of the hrst kind for the 
unknown kernel. Hence, with two measures, the identihcation of the kernel, 
which in principle is a nonlinear problem, is reduced to the solution of a 
mildly ill posed but linear problem. 

1 Introduction 

Let w = w{x,t), t > 0 and x G 12 (a bounded region with smooth boundary). 
The equation with persistent memory (A is the laplacian) 



( 1 . 1 ) 


is encountered in thermodynamics of systems with memory [w is the temper¬ 
ature) in nonhckian diffusion {w is the concentration) or in viscoelasticity. In 
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this case w(t) represents the displacement and the equation is usually written 
as 

w" = N{0)Aw{t) + [ M{t- s)Aw{s) ds, M{t) = N'{t). (1.2) 

Jo 

The vector ^ 

q{x,t) = / N{t — s)'Vw{x,s) ds 

Jo 

is the (density of the) flux at time t and position x in thermodynamics 
or nonfickian diffusion. In viscoelasticity, q'{x,t) is the traction. Let n{x) 
be the exterior normal to dfl at the point x. The products q{x,t) ■ n{x) 
or q'{x,t) ■ n{x), x G dfl, represent respectively the flux throughout the 
boundary or the traction on the boundary. In both the cases, they are 
observable quantities. 

Especially in viscoelasticity, N{t) is called the relaxation kernel. In this 
paper we assume that it is a smooth function (see below for details) and in 
the applications we mentioned above, it must be N{0) > 0. 

It is a fact that the kernel N{t) is a material property. Only its overall 
properties can be deduced from general principles (see for example M) and 
the specific function N{t) has to be identihed from experiments taken on 
the system. The “experiments” usually considered consist in the observation 
of the traction/flux on the boundary of a suitable sample of the material, 
subject to known excitations. The memory kernel N{t) is identified from 
such observations. 

Remark 1 The integrals in fll.ll) should extend from the time the material 
has been “created”, conceivably —oo. We assume that the measures are taken 
after a certain time to sufficiently far from the “creation” of the material and 
that the material was kept undisturbed for a long time, so that the effect of 
the memory of the previous history is negligible. This time to we rename 0. 

The shape of sample of the material usually does not influence the memory 
kernel otherwise, N would depend also on the position x (see i) so that we 
can try to determine N{t) from samples which have simple geometries. As 
in [1] and most of the references (commented in Sect. 11.21) . we shall use 
samples which are in the shape of a slab, so that we can assume that hi is a 
segment and we can normalize 

n = (0, vr). 
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The identification of N{t) is a nonlinear problem since if N{t) is not yet 
known, the function w{t) cannot be computed and the problem consists in 
the identification of the pair {N{t),w{x,t)). This pair enters in a nonlinear 
way in the equation. In spite of this, we shall see that using two measures 
we can construct a linear Volterra integral equation solely for N{t). So, the 
identification of N{t) is reduced to a linear problem. This is a first difference 
with existing literature. A second difference is as follows: in many papers, 
especially in engineering applications, the identihcation of N{t) is done as 
follows: it is assumed that the kernel belongs to a certain class of functions 
which depend on a finite number of parameters. The special kernel N{t) 
assumed to be the “true” kernel is the one which minimizes the discrepancy 
between a hnite number of experimental data and the corresponding theo¬ 
retically computed values. In contrast with this, we don’t assume any known 
a priori information on the functional form of N{t), a part smoothness and 
the condition A^(0) > 0. 

Due to different applications and terminology, for consistency and clarity 
we state that w will be called “temperature” and q will be called the “flux” 
(note that the usual minus sign in the definition of the flux or the traction 
has no influence in the following arguments and it is ignored). 

1.1 The assumption on N{t) and informal description 
of the procedure 

The hrst of our assumptions is that the kernel N{t) is of class (so that 
the memory kernel M{t) = N'{t) of the viscoelastic system is of class as 
often assumed in this kind of problems, see m)- 

It is seen from the form fll.2l) of Eq. fll.ll) that this equation is a per¬ 
turbation of the wave equation since we explicitly assume (? = A^(0) > 0. 
A consequence of A^(0) > 0 is that signals propagate in the body with finite 
velocity c = ^yN{0). So, whether this assumption holds is easily checked and 
the value of c is easily computed: just measure the velocity of propagation 
of the signals in the sample of the material. 

Our identification procedure applies in the case when c > 0. Solely to 
present simpler formulas, we assume that a first, simple, step is the identifi¬ 
cation of the value c = ^A^(O). We assume that this step has already been 
done. Then, the change of the time scale t ^ t/c transform N{t) to the new 
kernel Ni(t) = (l/c)A^(t/c). Note that A'i(O) = 1. So, we assume: 
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• the velocity c > 0 of signal propagation has been identihed. 

• the transformation of the time scale has been done already and we re¬ 
name as N(t) the function Ni(t) = {l/c)N(t/c). So, we study Eq. fll.ll) 
(i.e. (HID) with iV(0) = 1. 

Now we describe informally the identihcation procedure. We recall the 
goal: we fix any time T and we want to identify N{t) for t G [0,T] (we 
already know iV(0) = 1). 

1. We perform the following two measures, both on intervals of duration 
T (it is not restrictive to rename such intervals as [0,T]): 

(a) the temperature w{x,t) at the ends x = 0 and x = vr is kept 
equal zero, ta(0,t) = 0, w{7r,t) = 0. We impose a special initial 
condition tc(x, 0 ) = ^(x) and we read the flux y^{t) = g( 7 r, t) (note 
the lower case y). We use a suitable initial condition (a ramp) 
in order to identify an intermediate function K{t) which is then 
used in the identification of the memory kernel. 

(b) We impose null initial condition and w{7r,t) = 0 but now we 
excite the system from the left end: we impose w{0,t) = f(t). 
We measure q{7i,t). For clarity, the flux due to the boundary 
displacement / (with zero initial condition) is denoted Y(t) = 
Y^it) (note the upper case). 

In this second measure we don’t need to use any special function 
/(t), in particular f{t) is not related to the initial condition used in 
the hrst measure. But in practice the function f[t) will be smooth. 
For continuity we have to take /(O) = 0 and fit) “saturates” in 
a short time, i.e. it approaches lim 4 _j.+oo fit) (which in practice is 
finite) quite fast. 

2. We derive an algorithm for the identification of Nif) using these two 
observations. 

We note that the hrst step might be done by imposing w{x,t) = 0 and 
w'{x,0) = r]{x), with a special initial velocity r]{x) (integrating ( 11 . 21 ) we see 
that in this case Eq. (11.11) must have an additive term r]{x) on the right hand 
side). For dehniteness, we study rj = 0 and 7 ^ 0 but in practical applications 
to viscoelasticity = 0 y ^ 0 might be easier to impose. 
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Now we describe the organization of the paper. In Sect. 11.21 we discuss 
few references on the identihcation problem and in Sect. [2] we present the 
preliminary properties of system fll.ip we shall needed. The identification 
algorithm is in Sect. [3l 

It is possible to identify N{t) using a similar algorithm, which uses two 
measures with different initial conditions (and zero temperature at the bound¬ 
ary). We believe that this algorithm is not practical since to impose initial 
condition is more difficult then to impose a boundary temperature (or de¬ 
formation). However, for completeness, this second algorithm is described in 
Sect. 13.31 

1.2 Comments on previous references 

Due to the importance of this problem, reference on kernel identihcation is 
enormous and both Engineers and Mathematicians addressed the problem, 
along different and scarcely related lines. Most of the available results can 
be described as follows. An approach to kernel identihcation, seemingly hrst 
proposed in [13] and then followed by a wealth of related papers (see also the 
books mm for references) is as follows: it is noted that the identihcation of 
N{t) on the basis of additional measures is in fact the identihcation of the pair 
{w(t), N(t)) and w(t) solves an integrodiherential equation. Now the trick 
is to compute derivatives of the measured output (of course when the data 
are sufficiently smooth) and to derive a second integrodiherential equation 
for N{t). In this way we get a nonlinear integrodiherential equation for the 
pair N{t)) and, in general using hxed point theorems, it is possible to 

prove existence of solutions of this systems. 

This leads to the identihcation in particular of N{t) (according to the 
assumptions, either locally, for a short time, or globally, as in [T3]). 

The identihcation requires only one measure, taken assuming known (and 
smooth) initial condition (or affine term). 

The approach that is most common in engineering literature is as follows. 
It is assumed that the kernel N(t) belongs to a certain class of functions, 
identihed by a hnite number of parameters. The classes that are most often 
considered are Dirichlet sums, i.e. 


K 



(1.3) 


5 






as for example in [H [SI |2H 125] (here a > 0, 6 > 0) or a combination of Abel 
kernels 

= ( 1 - 4 ) 

k=0 ^ 

sometimes a combination of both (see for example mm)- 

The class of the kernels (ll.Sp or (ll.4|l have a “physical” sense since they 
are encountered when solving ordinary differential equations of integer or 
fractional order and these classes of kernels are reported to £t sufficiently 
well experimental data, but not in every application, see [23] . 

In order to identify the parameters, a sample of the material, usually in 
the form of a slab, is excited using a known initial or boundary deformation. 
The stress relaxation is measured on (a part of the) boundary and the pa¬ 
rameters (ofe, bk, sometime also the exponents (3k) are chosen so to fit better 
the experimental measures (usually minimizing a quadratic index). This can 
be done using a single measure. Multiple measures can be used to improve 
accuracy or to remove spurious frequencies. 

Different algorithms works in the frequency domain, using Fourier or 
Laplace transformations, as for example in na. 

The use of a single measure is common to the algorithm presented in [T3l 
m and related papers. Here instead we assume that two different measure 
are performed and, from these measures, the kernel N(t) is reconstructed. We 
don’t assume that the kernel belongs to any specific class of functions, a part 
regularity and the physical property A^(0) > 0. In fact, no further properties 
of the kernel are in principle needed for the reconstruction algorithms (not 
even positivity for every time!) in spite of the fact that such conditions 
are imposed by thermodynamics and can be used to improve the results, 
see Sect. 13.41 Note however that we are assuming that N{t) is smooth for 
t > 0 and this rules out Abel kernels. The extension of our algorithm to Abel 
kernels will be pursued in future researches. 

2 Preliminary properties of system ( ll.ll) 

We present a representation of the solutions of system fll.ip when D = (0, tt), 
taken from [niiisiiis] (see [2I1I22] for more general results, in particular for 
extensions to the case dimfl > !)• 

We consider the sequence of the functions 
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'n - \/ -smnx, 


n e N. 


It is known form the theory of Fourier series that {0n} is an orthonormal 
basis of n) and that 

A(pn{x) = = -n^4>n{x) . 

We expand the initial condition ^ and the solution in sine series: 

PTT PTT 

in = I sin nx dx , Wn{t) = / w{x,t)smnxdx 


+ 00 


= ^^nsmnx, 


n=l 


+00 


w{x, t) = Wn{t) sin 


nx 


n=l 


w'^ = -n^ N{t-s)wn{s)ds + n N{t - s)f{s) ds, Wn{0)=U- 


>0 


>0 


This is a scalar integrodifferential equation whose solution can be represented 
as 

Wn{t) = U^nit) + n Zn{t-s) N{s - v) f (v) dv ds 

. Jo Jo 

where Zn(t) solves 

^n{t) = f - s)Zn{s) ds , Zn{0) = 1 • (2.1) 

Jo 

The transformation (^,/) ^ w is linear and continuous from L^(0,7r) x 
L2(0,T) to C'([0,T];L2(0,7r)) x C^i[0,T]; H-\0,7r)). 

We need the following observation; 

Lemma 2 Let T > 0. There exists M = Mt such that for every t G [0,T] 
and every n we have: 


4(^) 


n 


< 


M 


n 


a = --iV'(O), 


( 2 . 2 ) 
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Proof . This inequality follows for example from the equality (5.4) in [T7] (see 
also MDE]). We must be careful with the notation. With a = —N'{0)/2 
let N{t) = e^"*iV(t) so that iV(0) = 1 and iV'(O) = 0. The functions denoted 
Znit) and N{t) in [T3, formula (5.4)] are respectively and N{t) = 

(see [13 formula (5.1)]). So, we have the following equality for the 
function z„(t) in fl2.ll) (here (3n = — a^): 


Zn{t) = e "‘cos/3nt + ^6 siu ^ / e ‘^°'^Zn{s) N'{t - s) + 


2 rt 


A 


A. 


n ^ 0 


(*t—S 


+ 


aN'{t — s — r) — N"{t — s — r) 


e"'’ cos dr 


ds. 


We compute the derivatives of both the sides, using A^'(O) = 0- We get the 
equality 


cy 

z'nit) + Ane"“^ sin PJ, = sin Pnt- 

Pn 


/ e ^"^z„(s) <! iV"(t - s) - iV"(0)e"^* cos Pn{t - s) 
Pn Jo '' 


2 pt 


ft — S 


+ 


aN"{t — s — r) — N'"{t — s — r) e“^ cos PnV dr > ds . 


The result follows by dividing both the sides with n and nothing that if 
t G [0,T] then there exists M such that 


I sin Pnt — sin t\ < — . i 
n 

We are going to use these properties of the system for the identihcation 
of the kernel N{t), assumed unknown. Note that these same properties have 
also been used to solve a problem of source identihcation, see [siEni. 


3 The identification procedure 

Now we give the details of the two steps of the identihcation procedure de¬ 
scribed in Sect. o 













3.1 A first measure of the flux, using initial tempera¬ 
ture 

Using a special initial temperature, we can identify an intermediate function 
K(t) which will be used in the second step of kernel identihcation. 

Let w{x, 0) = ^ ^„sinnx G -L^(0, vr) so that G P and 


-hC50 

'^w{x,t) = 'Y^inZn{t) sinnx G U([0,T];L^(0,7r)) , 

n=l 
+ 00 


'^w^{x,t) = 'Y^{n^n)zn{t)cosnx G U([0,T];iL ^(0,7r)). 





72=1 




So we have 







TT 

= 

-\-oo 


rt 


+00 ^ 


E«” 

n 

/ N{t - s)Zn 

'0 

(s) ds 

cosnx = - 

cos nx 


72=1 



72=1 

(3.1) 

and 










TT 

= = 

-l-CXD 

1 




(t) 

-E 

n 

(3.2) 


Using the inequality fl2.2p . we see that y^{t) G L^(0,T) for every T > 0 and 
every ^ G L^(0, vr). If in particular we choose 

1 1 

^ 0 ( 2 ^) =-(tt — x) = > — sinnx (3.3) 

2 n 

72=1 

then the series in fl3.2p is 

+00 .. 

72=1 

The function we wanted to identify is 

+00 .. 

A'(*) = ■ P,4) 

72 = 1 

This is the hrst measure we need in the identihcation process. 
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3.2 Measure of the flux due to a boundary excitation 

Now we combine the identification of K{t) obtained in the previous section 
with an observation of the flux at x = tt due to a nonzero temperature applied 
to the left hand x = 0 (and zero initial condition). Hence we consider the 
following system: 


w'{t) = J N{t — s)Aw{s) ds , 0) = 0, 

w{0,t) = f{t), w{7r,t) = 0. 


(3.5) 


We choose any input / such that /(O) = 0, (this is consistent with the 
fact that the initial condition is zero). So we have 

fit) = [ ais) ds 

Jo 

and we explicitly assume 

^?( 0 ) = /'( 0 )^ 0 . 


We use 


71 /“* r 

—w{x,t) = 2 ^{smnx) / /(s) n / N{r)zn{t — s — r) dr 

2 Jo L Jo 


n=l 
+ 00 


ds = 


y^fsinnx) — [ 

I — (sin nx) 


1 /■* d 

fis)—Znit - s) ds 
n In ds 


\n=l 


rt pt 

g{r) dr- gis)znit - s) ds 


(3.6) 


The function {x,t) i—)■ w{x,t) belongs to C'([0, T]; L^(0, vr)) so that (x,f) 
Wxix,t) belongs to ^([O, T]; ih“^(f2)), 


71 

2 


Wxix,t) 


pi \ - 1-00 pi 

/ (^(r) dr I cosnx — (cosux) / g{s)znit — s) ds = 

'to J n=l n=l 'to 

\ H-oo pi 

g{r) dr ) - ^ (cosnx) / g{s)zn{t 

J n=l Jo 


7r6{x) ( / g{r) dr ]- 


— s) ds. 
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Here 6{x) is the Dirac’s delta (supported at x = 0) and so when observing 
the flux at x = TT we get 


7r ^/\ TT, ^ 1 

^YHt) = 


N{t — s) J g{r) dr dsj — 

pt pt — S 

9 {s) N{t - s - r)zn{r) dr ds = 

n=i Jo Jo 

(/ lo ^ Jo 


N{t — s) / g{r) dr ds I + / g{t — s) 


H-cxd 


n=l 




= / »(< - s) 1 K{s) - ^ f N{r) dr ) ds 


ds 

(3.7) 


where, from 03.41) . 


A-(s) = 


+00 




n=l 


We have seen in the previous section that K{t) can be estimated from the 
measure of one boundary flux. So, at this point, K{t) is a known function. 
Formula 03.7p explains its significance: K{t) — 1/2 N{s) ds is the impulse 
response of the system with boundary input / at x = 0 and the flux observed 
at X = TT. 


Remark 3 In fact, the 5-function is an artifact, due to the fact that for¬ 
mula 03.61) gives an odd extension of the solution which is periodic of period 
27r and so the formula introduces a jump at f = 0 (because f(t) ^ 0) and 
this produces a 5 function at x = 0 and also at every multiple of ‘2'k. This 
fact has been discussed in lais]. I 


The function iV [t) can be identified from formula 03.7p which is a Volterra 
integral equation of the first kind for 

M{t)= f N{s)ds. 

Jo 
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Although this may not be the best way from the point of view of numerical 
computation, we may note that if (? G hence / G then Y(t) is 
differentiable and upon computing the derivatives of both the sides of fl3.7p 
we get: 


g{0)M{t)+ [ g'{t-s)M{s) ds = 

Jo 

= ‘^g{0)K{t) + 2 ^ g\t - s)K{s) ds - TT • 

We have chosen an input signal / such that /'(O) = 5 ^( 0 ) 7 ^ 0. So, this is a 
Volterra integral equation of the second kind (but note that the right hand 
side requires the numerical computation of Y'{t)). 

Once this Volterra integral equation has been solved, a last step of nu¬ 
merical differentiation identifies N{t). 

Remark 4 We note that 

• If 5 '( 0 ) = 1 and g'{t) = 0 then we get an explicit formula for M{t) 
but in this case f{t) = t. A steadily increasing temperature cannot be 
applied. In practice, f(t) has to saturate, as for example f(t) = tj 

or 1 — e“*. 

• In practice it does not seems difficult to impose a variable temperature 
(or traction) on the boundary. It may be more difficult to impose a 
special initial condition but we note that the initial condition we need 
in the identification of A' is a ramp, which may be achieved by letting 
the equilibrium be reached between two different temperatures of the 
ends of the sample. 1 

3.3 A different identification procedure 

We present now a second identification algorithm which uses two different 
measures, obtained using two different, but related, initial conditions. In the 
algorithm we proposed above, the most questionable step seems the use of 
a special initial condition. It is even more difficult to impose also a second, 
related, initial temperature. So, the procedure we present now may not have 
a practical interest, but it seems mathematically interesting. 
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We impose an initial temperature ^(x) = (13■2p : 

AU,. ,w4W 


»'(() = -5Z(-i)"«" 


n=l 


n 


Let us compute the convolution N *y^ (we recall that the kernel N[t) is still 
unidentihed): 


+ 00 


(-ir 


N{t - s)y^{s) ds = - V- —^n / N{s)z'^{t - s) ds . 

“ n In 


(3.8) 


n=l 


This formula is easily interpreted: let us consider the second initial condition 


H-OO 


= Y1 


sm nx = X 


n=l 






r) dr ds 


and let us compute y'^{t) from the intermediate equality in (13.111 (where we 
put X = TT and we replace with We get 

+ 0O ^ 

/ N{s)zn{t - s) ds 


so that 




N{s)z'^{t — s) ds (3.9) 


\n= 
/ +cx> 


(Z j - I P 

The functions y^ and y'^ can be measured and so we hnd a Volterra integral 
equation for the unknown function N{t). The Volterra integral equation is 
of the second kind, i.e. its solution is a well posed problem, if ^ is such that 


+ 00 

E 

72=1 




n 


^0. 


In spite of this, the identihcation problem is ill posed since the procedure 
requires the numerical computation of the derivative of y^{t). 
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3.4 A toy simulation 

Whether the identification procedure we presented is really useful is a matter 
to be settled by the Engineers on the basis of concrete experiments. Here we 
present a toy simulation on the basis of theoretical computations. We assume 
that N{t) is a linear combination of exponentials, as it is often assumed in 
applications, and for simplicity we confine ourselves to the case that N{t) is 
the combination of three exponentials, 

N{t) = (3.11) 

where the coefficients and the exponents are positive. Then, the equa¬ 
tion of Znit) is 

^hit) = -1- Zn{s) ds . 

Jo 

It is easy to reduce the computation of Zn{t) to the solution of an ordinary 
differential equation. In fact, Zn{t) is the first component of the solution 
{z,Wi, 102 , 103 ) of the following problem: 


( z' = -n^{ciWi C2W2 + C3W3) 

( m \ 


( ^ \ 

1 W[ = -TiWi + Z 

Wi(0) 


0 

1 W2 = -T2IO2 + ^ 

m;2(0) 


0 

[ Wg = -T3W3 + Z 

\ w^3(0) / 


VO/ 


So, the functions Zn(t) can be easily computed and we can compute also 
Y^{t) and K{t). Then we choose f{t) = 1 — and we reconstruct the 
kernels M{t) = J^N{s) ds and N{t). 

Before we present the simulation we note that: 

• the reconstruction of the relaxation kernel is in two steps: first we 
identify M{t) and then a next step of numerical differentiation gives 
N{t). The identification of M{t) is a deconvolution problem which can 
be solved, essentially with the same results, with Tikonov regularization 
(as described in m) or Lavrientev method (see 0171). Both these 
methods introduce fast oscillations but it is known in this problem that 
the oscillations are an artifact due to the amplification of the noise. In 
fact, M{t) is increasing (since for physical reasons N{t) > 0) and, in 
the physical case that N(t) is decreasing, M(t) is concave. Hence these 
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oscillations can be attennated for example by compnting the average 
of the reconstruction of M{t) on nearby steps. This has been done in 
the simulation below. 

• in the simulation we present below we choose to reconstruct M{t) using 
Lavrentev method (followed by an average, to remove fast oscillations). 
Instead we use Tikonov regularization in the next step of the recon¬ 
struction of N{t) from M[t) since Lavrentev method is more efficient 
when it is known that the function to be reconstructed is 0 at the initial 
time (and since Lavrentev method is more sensitive to disturbances). 

• Also the function N{t) can be averaged to reduce the effect of errors 
since in practice we know that it is a decreasing (and often a convex) 
function. Taking into account the known value iV(0) this last step (not 
shown in the plots) will also improve the reconstruction near t = 0. 

The special case we test is as follows: the time interval is [0, 5] and we 
assume 10 measures every time unit. We choose 

N{t) = (l/10)e-*/2 + (i/5)e-2* + (1/2)6-^' (3.12) 

(i.e. we normalize so to have iV(0) = 1). The functions K[t) and Y^[t) are 
approximated (and corrupted by an error of 1% at every “measure”). Then, 
we “reconstruct” both M{t) and N{t). Figure [U (left) presents the plots 
of M{t) and its reconstruction from Lavrientev method (the plot where the 
amplification of the disturbances is evident) and its average (at every step 
tree values on the left and on the right have been averaged). The averaged 
reconstruction of M{t) is fed to the Tikonov algorithm in order to compute 
the numerical derivative of M{t), i.e. N{t), which is plotted on the right. 
The penalization parameters of the Lavrientev and Tikonov algoritms are 
respectively 0.1 (equal to the interval of time between consecutive measures) 
and 0.01. 
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Figure 1: Left: the kernel M{t) (the function, the reconstruction and its 
average); right: the reconstruction of the relaxation kernel N{t) (and its plot 
from the expression fl3.12p 
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